function cfdTimeAdvanceRK4
%--------------------------------------------------------------------------
%  Written by LiuHaHa @ NWPU, 2022 01 
%  Contact me at: liuzhikan@mail.nwpu.edu.cn
%==========================================================================
% Routine Description:
%   This function advance time using 4th order Runge-Kutta
%--------------------------------------------------------------------------
global cfdFlow slvpara;
theNumberOfElements = cfdFlow.mesh.numberOfElements;
phyTemp = cfdFlow.Phy;
phyPre = cfdFlow.Phy;
elementVolumes = cfdFlow.mesh.elementVolumes;
cfdFlow.phyPreOutput = cfdFlow.Phy;

if (slvpara.flx_order == 1)
    RKa(1)  = 0.0833;      
    RKa(2)  = 0.2069;      
    RKa(3)  = 0.4265;     
    RKa(4)  = 1.0000;      
else
    RKa(1)  = 0.1084;
    RKa(2)  = 0.2602;
    RKa(3)  = 0.5052;
    RKa(4)  = 1.0000;
end

for subIrt=1:4
    
    cfdCalculateFluxAndEigenvalue;
    %
    cfdCalculateResidual;
    
    if (subIrt == 1)
        cfdCalculateDeltaT;
    end
    
    dt = cfdFlow.dt;
    resd = cfdFlow.Resd;
    for idim=1:slvpara.ND+2 
        phyTemp(1:theNumberOfElements,idim) = phyPre(1:theNumberOfElements,idim)...
            - RKa(subIrt)*dt*resd(1:theNumberOfElements,idim)./elementVolumes(1:theNumberOfElements);
    end

    cfdFlow.Phy = phyTemp;
    
end

cfdFlow.PhyTotalTime = cfdFlow.PhyTotalTime + cfdFlow.dt;
















